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ABSTRACT 


A computer  program  developed  specifically  for  computing  the  heat  loss  from 
a direct  buried  underground  heat  distribution  system  based  on  the  measured 
values  of  soil  thermal  conductivity  and  the  earth  temperatures  over  the 
underground  pipes  is  presented.  The  heat  loss  rates  and  locations  of  the 
underground  pipes  were  calculated  by  statistically  determining  the 
parameters  in  a theoretical  expression  for  the  underground  temperature 
field  around  a two-pipe  system  using  the  nonlinear  least  squares  method. 
Sample  calculations  based  on  two  sets  of  test  data  obtained  from  the  scale 
model  experiments  are  presented,  and  the  results  obtained  from  this 
computer  code  implemented  on  a microcomputer  are  compared  with  those  by  the 
DATAPLOT  software  package  installed  on  a mainframe  computer. 


Key  Words:  Computer  program,  district  heat  distribution  systems,  nonlinear 
least  squares  fitting,  underground  heat  distribution  systems. 
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Introduction 


The  heat  loss  from  underground  heat  distribution  systems  accounts  for  a 
considerable  amount  of  the  fuel  energy  consumed  to  transport  hot  water  or 
steam  from  a power  plant  to  various  buildings.  A recent  field  survey 
conducted  by  Pan  Am  [1]  on  Department  of  Defense  and  Veteran  Administration 
installations  clearly  demonstrated  that  the  thermal  performance  of  the 
direct  buried  underground  system  is  degrading  due  to  deteriorating  pipe 
insulation,  and  corroding  pipes  and  conduits  resulting  in  ground  water 
penetration.  The  field  survey  also  indicated  that  the  repair  costs  for 
direct  buried  conduit  systems  are  higher  than  for  concrete  shallow  trench 
systems  due  to  extensive  work  required  to  excavate  and  cut  open  the  conduit 
pipes  to  locate  and  repair  failed  line  segments.  Accurate  and  simple  in- 
situ  methods  for  determining  the  heat  loss  from  the  underground  sytem  are 
needed  to  provide  information  about  the  existing  thermal  performance  of  the 
insulation.  The  heat  loss  from  the  system  under  actual  use  conditions  is 
crucial  for  determining  if  the  system  requires  either  repairing,  or 
replacing  its  defective  segments. 

The  thermal  probe  method  [2]  developed  at  NBS  is  a reliable  and  easy-to-use 
instrumentation  system  for  evaluating  the  performance  of  direct  buried 
underground  systems.  This  non-destructive  measurement  technique  is  based  on 
the  heat  transfer  theory  relating  the  earth  temperature  profile  over  the 
underground  pipes  to  the  heat  loss  values,  provided  that  the  thermal 
conductivity  of  the  surrounding  soil  and  the  separation  distance  between 
the  pipes  are  known.  The  amount  of  heat  loss  from  the  underground  pipes 
and  their  locations  were  derived  statistically  by  determining  the 
parameters  in  a theoretical  expression  obtained  from  a mathematical  model 
using  the  nonlinear  least  squares  technique  along  with  the  gathered  data. 
A large  software  package  DATAPLOT  [3]  developed  for  statistical  analysis  is 
available  for  linear  and  nonlinear  fitting,  graphics  and  data  correlation. 
This  package  needs  memory  storage  requirements  of  at  least  3 Mbytes  and  has 
been  implemented  on  mainframe  computers.  For  nonlinear  least  squares 
fitting,  this  software  package  makes  use  of  a Levenberg-Marquardt-Morrison 
algorithm  [4-6]  with  finite  difference  approximation  for  the  first  partial 
derivatives  of  the  non-linear  function  with  respect  to  the  parameters  to  be 
determined . 

Microcomputers  are  widely  used  as  a tool  to  solve  complex  mathematical  and 
logical  problems  in  scientific  and  engineering  applications.  With 
prodigious  techno Igical  growth  and  advances,  microcomputers  are  far  easier 
to  operate  and  cheaper  to  upgrade  the  systems  for  computing  speed,  power 
and  versatility  than  the  large  main  frame  and  mini-computers.  The  popular 
trend  of  using  micro-computers  in  experimental  work,  such  as  automatic 
control  of  instrumentation  systems  and  data  acquisition  and  processing,  is 
anticipated  to  continue  since  micro-computers  carry  out  various  functions 
at  low  cost. 

This  report  is  a summary  of  work  in  progress  and  desribes  a computer 
program  developed  specifically  for  calculating  the  heat  loss  from  an 
underground  heat  distribution  system  based  on  the  earth  temperature  and 
thermal  conductivity  data  using  the  non-linear  least  squares  method.  A 
two-pipe  configuration  for  the  underground  system  is  selected  for  numerical 
treatment  since  most  of  heat  distribution  sytems  consist  of  heat  supply  and 
return  pipes  of  different  sizes  and  pipe  temperatures.  The  computer 
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program  utilizes  the  earth  temperature  data  points  to  solve  for  parameters 
in  a steady  state  heat  conduction  equation  for  the  underground  temperature 
field  around  a two-pipe  system.  The  predicted  values  of  these  parameters 
yield  the  calculated  values  of  the  heat  loss  rates  and  locations  of  the 
pipes  as  the  outputs  of  the  computer  program.  In  the  field,  the  locations 
of  the  underground  pipes  may  not  be  known.  The  outputs  from  this  computer 
code  can  also  be  used  to  provide  a quick  means  of  evaluating  the 
experimental  data  on  soil  temperature  and  thermal  conductivity  during  field 
measurements.  The  algorithm  is  basically  the  Levenberg-Marquardt-Morrison 
method  [4-6]  for  the  fitting  of  the  multivariable,  non-linear  earth 
temperature  equation  by  unconstrained,  unweighted  nonlinear  least  squares 
with  the  first  partial  derivatives  derived  analytically. 

Also  presented  in  this  report  are  sample  calculations  based  on  two  sets  of 
test  data  obtained  from  the  scale  model  experiments.  The  simulated  scale 
model  of  an  underground  system  was  essentially  an  insulated  five  sided  box 
consisting  of  two  electrically  heated  pipes  buried  in  dry  sand.  The 
electrical  power  inputs  to  the  pipes  and  the  sand  temperatures  in  the 
vicinity  of  the  buried  pipes  were  measured  during  the  tests. 


2.  The  Earth  Temperatures 


The  underground  temperature  field  caused  by  two  buried  pipes  under  steady- 
state  conditions  can  be  expressed  as  follows  [2]: 


T 


2 Qi  (X-b.O^+Cr+ai)^ 

Z — In  [ ^:]  + T„ 

i=l  (X-bi)^+(Y-ai)^ 


(1) 


where  T = the  temperature  of  the  soil  at  a given  location 

= the  strength  per  unit  length  of  the  i-th  pipe  (heat  source  or 
sink) 

k = the  thermal  conductivity  of  the  soil 
X,Y  = the  position  of  any  arbitrary  point  in  the  temperature  field 
^i»^i  " horizontal  distance  and  vertical  depth  of  the  center  of  the 

i-th  pipe 

Tq  = the  soil  temperature  without  the  pipe's  present 


This  non-linear  multivariable  function  can  be  solved  to  give  the  pipes' 
heat  losses  (Q]^,Q2)j  locations  (bpb^)  and  depths  (apa2)  using  the  method 
of  non-linear  least  squares,  provided  that  the  earth  temperature  and 
thermal  conductivity  data  are  available.  In  order  to  improve  the  estimate 
of  pipe  heat  loss,  one  of  the  unknown  parameters  is  removed  with 
introduction  of  a known  separation  distance  between  the  centers  of  the 
pipes.  This  separation  distance  can  be  obtained  from  either  the  pipeline 
layouts  in  architectural  drawings  or  actual  measurement  with  the 
underground  pipes  located  in  the  nearest  accessible  manholes.  This 
function  derived  by  the  heat  transfer  theory  has  five  unknown  parameters 
which  must  be  adjusted  until  an  optimum  fit  is  obtained  between  the 
predicted  temperature  and  the  measured  temperature  data. 


3.  Algorithm  for  Non-Linear  Least  Squares 

The  algorithm  for  non-linear  optimization  of  the  multivariate  functions  is 
basically  the  Levenberg-Marquardt-Morrison  method  [3-5].  The  function  is  a 
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sum  of  squared  terms,  which  are  a function  of  one  or  more  parameters,  to  be 
minimized : 


N 

S(v)  = Z [Ri(Vi,V2, 


i=l 


• • • , 


(2) 


and  the  residual  of  the  predicted  and  observed  values  is 


R^Cv)  = F(V,Z)  - Y 


(3) 


where  V = vector  of  the  unknown  parameters 
F = a non-linear  or  linear  function 
Z = vector  of  one  or  more  independent  variables 
Y = vector  of  observed  values  or  a dependent  variable 
N = total  number  of  observed  values 
m = the  number  of  parameters  to  be  determined 

The  vector  of  parameter  estimates  after  r+1  iterations  can  be  obtained  from 


y(r+l)  . v(r)  + IJ 


(4) 


By  solving  the  following  equation  derived  from  the  Taylor  series  expansion 
of  a function  [3-5],  the  step  length,  h,  to  minimize  the  non-linear 
function  can  be  obtained 


(5) 


(A'A  + p B)h  = -A'R 


where  A = the  matrix  of  first  partial  derivatives  of  R^'s  with  respect  to 


A'  = tne  transpose  of  matrix  A 


p = the  Levenberg  parameter 

B -a  diagonal  matrix  with  the  same  diagonal  elements  as  A'A 
R = vector  of  R^'s 

In  this  computer  program  for  non-linear  regression,  all  of  the  first 
partial  derivatives  are  derived  analytically  instead  of  by  finite 
difference  approximation  as  used  in  the  DATAPLOT  software  package  [3]. 

The  set  of  simultaneous  equations  derived  from  equation  5 is  solved  for  h 
by  using  Householder  transforms  [7]  to  reduce  the  matrix  A to  upper- 
triangular  form.  The  reduced  upper-trangular  matrix  is  then  added  to  the 
square  root  of  B. 

In  any  iteration,  the  Levenberg  parameter  is  multiplied  by  the  factor  DECR 
if  the  sum  of  squares  of  the  residuals  is  sufficiently  reduced  at  the  first 
attempt,  and  by  the  factor  EXPND  if  not.  A new  step  length  h is  then 
calculated . 

A sufficient  reduction  in  the  sum  of  squares  is  defined  as  the  condition 
when 


PSI  = (0.5)(S§  - sf)/(sg  - S^)  < 10“*^ 


(6) 
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where  Sq,  S2,  and  S2  square  roots  of  the  sum  of  squares  at  the 

start  of  the  iteration,  at  the  finish  after  taking  the  step  h,  and  the 
residual  sum  of  squares  at  the  finish  of  an  iteration,  respectively. 

The  test  for  convergence  of  the  sum  of  squares  of  the  residuals  is 
considered  to  be  satisfied  when 

(So-S2)/(l+So)  < TOL 

where  TOL  is  a tolerance  level  specified  in  testing  for  convergence. 

4.  Program  Description 

The  computer  code  consists  of  a main  program  (UHDS)  and  two  subroutines 
(LMMNLF  and  FUNVAL).  The  main  program  handles  all  input  data.  Output  is 
from  the  main  program  and  LMMNLF.  Subroutine  LMMNLF  is  based  on  the 
Levenberg,  Marquardt  and  Morrison  algorithm  [2]  to  find  a minimum  of  a sum 
of  squares  which  is  a function  of  one  or  more  parameters.  This  subroutine 
has  been  modified  to  be  applicable  to  a nonlinear  function  of  more  than  one 
independent  variables.  This  subprogram  called  from  the  main  program 
performs  and  coordinates  all  calculations,  and  provides  the  output  of  final 
parameter  estimates.  Subroutine  FUNVAL  is  used  to  evaluate  the  nonlinear 
function  and  to  carry  out  calculations  of  its  first  partial  derivatives.  A 
listing  of  this  computer  program  is  given  in  Appendix  B. 

(1)  Description  of  Variables: 

X Vector  of  the  parameters  in  the  nonlinear  function 

F Vector  of  function  values  (residuals)  whose  sum  of  squares 

is  to  be  minimized 

A A matrix  containing  the  first  partial  derivatives  of  the 

function  with  respect  to  each  of  the  parameters 

SUMSQ  The  final  residual  sum  of  squares 

N Number  of  terms  in  the  sum  of  squares  or  the  number  of  data 

points 

NP  Number  of  the  parameters  to  be  determined 

TOL  Tolerance  for  residual  sum  of  squares 

EXPND  A factor  by  which  the  Levenberg  parameter  EPS  is  increased 
when  difficulty  is  experienced  in  reducing  the  sum  of 
squares 

DECK  A factor  by  which  EPS  is  decreased  if  the  sum  of  squares  is 

reduced  at  the  first  attempt  in  an  iteration 

ITS  Input:  Maximum  number  of  iterations  to  be  allowed 

Output:  Actual  number  of  iterations  performed 
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lER 

Input:  =0  No  print-out 

= 1 Print  parameter  estimates,  EPS,  PSI,  residual 
standard  deviation  after  each  iteration 

Output:  = 1 Successful  termination 

= 2 Maximum  iterations  exceeded 
= 3 EPS  exceeds  10^ 

= 4 Attainable  accuracy  has  been  reached,  TOL  was 
set  too  small 

IFL 

A parameter,  IFL  = 1 First  partial  derivatives,  RESID,  SUMSQ 

are  to  be  calculated 

IFL  = 2 Calculate  SUMSQ  only 

EPS 

The  Levenberg  parameter 

NIV 

Number  of  independent  variables 

NRDF 

Number  of  residual  degress  of  freedom 

SDRES 

Residual  standard  deviation 

PRED 

The  fitted  value  of  the  dependent  variable 

PSI 

A ratio  as  defined  in  equation  6,  to  be  tested  if  a 
sufficient  reduction  in  the  sum  of  squares  has  been  reached 

RESID 

The  difference  between  observed  and  fitted  values  of  the 
dependent  variable 

AK 

Soil  thermal  conductivity,  (Btu/h-f t-°F) 

DS 

Separation  distance  between  the  centers  of  the  pipes,  in 
inches 

NSETS 

Number  of  the  data  sets 

YY 

The  dependent  variable  (the  earth  temperatures)  (deg.  F) 

XX(I,J) 

The  independent  variables,  I = Data  point  number, 

J = 1 Horizontal  distance,  (inches) 

= 2 Vertical  depth,  (inches) 

= 3 Undisturbed  earth  temperature,  (deg.  F) 

= 4 7.958  X 10"2/ak 

=5  Separation  distance  between  the  pipe  centers, 
(inches ) 
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(2)  Input  Description 


A file  called  DATAFL  to  serve  as  the  input  data  file  is  needed  to  be  set  up 
prior  to  execution  of  the  computer  program.  This  file  contains  data  to 
specify  the  number  of  data  sets,  the  number  of  data  points,  the  soil 
thermal  conductivity  and  the  distance  between  the  centers  of  the  pipes  for 
a given  data  set,  the  measuring  locations,  the  observed  earth  and  the 
undisturbed  soil  temperatures,  and  the  initial  estimates  of  the  parameters 


to  be  determined.  Complete  descriptions  of 

specification  in  the  input  file  are  given  below: 

input  formats  for  each 

Record  No.  Format 

Contents 

1 (15) 

NSETS 

2 (15,  2F10.5) 

N,  AX,  DS 

3 (F5.1,  2(1X,  F4.1),  IX,  F5.1) 

YY(I),  (XX(I,J),  J = 
1,  3 

(Use  N records  of  this  form) 


N+3  (F10.5) 

Initial  estimate  of 
X(l):  Heat  loss 
rate  of  pipe  No.  1 

N+4  (F10.5) 

X(2):  Horizontal 
distance  of  pipe  No. 
1 

N+5  (F10.5) 

X(3):  Vertical 
depth  of  pipe  No.  1 

N+6  (F10.5) 

X(4):  Heat  loss 
rate  of  pipe  No.  2 

N+7  (F10.5) 

X(5):  Vertical 
depth  of  pipe  No.  2 

(Repeat  record  numbers  2 to  (N+7)  for  the  second  data  set) 
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5.  Sample  Calculations 


Two  sets  of  experimental  data  obtained  from  scale  model  experiments  were 
used  to  test  the  computer  code.  The  calculations  were  performed  on  an  IBM 
personal  computer.  The  computer  outputs  for  these  sample  calculations  are 
contained  in  Appendix  A. 

In  order  to  check  the  validity  of  the  calculation  scheme,  the  results 
obtained  from  this  computer  code  were  compared  with  those  by  the  DATAPLOT 
software  package  implemented  on  an  UNIVAC  1100/80  computer.  Table  1 
presents  these  compared  results  of  the  pipe  heat  loss  rates,  and  the 
horizontal  distance  and  vertical  depth  of  the  centers  of  the  electrically 
heated  pipes  for  both  sets  of  the  test  data.  In  general,  there  is  a small 
discrepancy  between  these  two  computation  schemes  on  the  calculated  values 
except  for  the  heat  loss  rate  of  pipe  No.  1 in  data  set  No.  1.  However, 
the  total  heat  loss  from  both  pipes  calculated  from  this  computer  code 
agrees  reasonably  well  with  that  obtained  from  the  DATAPLOT  software 
package.  It  can  be  noted  that  the  standard  deviations  of  the  heat  loss 
from  pipe  No.  1 for  both  data  sets  are  significantly  larger  as  illustrated 
in  the  computer  outputs.  This  may  be  due  to  the  variability  of  sand 
thermal  properties,  which  resulted  in  an  irregular  temperature  distribution 
in  local  areas  observed  experimentally,  especially  in  the  vicinity  of  the 
heated  pipes,  and  some  experimental  errors  of  the  sand  temperature  and 
thermal  conductivity  measurements.  In  Table  1,  the  corresponding  measured 
values  for  these  scale  model  experiments  are  also  listed  for  comparison. 
It  can  be  seen  that  the  calculated  values  are  generally  in  good  agreement 
with  the  measured  results. 

6.  Conclusions 

A computer  program  has  been  developed  specifically  to  carry  out 
calculations  for  estimating  the  heat  loss  and  the  locations  of  a pair  of 
underground  pipes  based  on  the  non-linear  least  squares  method.  This 
program  can  be  implemented  on  a micrcomputer  and  gives  proper  computing 
speed  and  adequate  accuracy  on  the  calculated  results.  The  soil 
temperature  and  thermal  conductivity  data  obtained  from  the  thermal  probe 
method,  and  the  separation  distance  between  the  pipes  are  required  as  the 
input  data.  The  calculated  pipe  heat  loss  rates  can  be  used  to  identify 
any  anomalies  which  may  be  indicative  of  defective  segments  of  a direct 
buried  conduit  underground  system.  The  outputs  from  the  program,  the 
residuals  between  the  measured  and  predicted  soil  temperatures  for 
instance,  can  serve  as  a rapid  means  for  evaluating  the  results  of  various 
measurements  on  the  test  site.  The  computational  scheme  and  the  input  data 
file  required  for  executing  this  computer  code  are  described,  and  the 
computer  outputs  for  two  sample  cases  are  presented. 
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Table  1 

A Comparison  of  Calculated  Results  from  this  Computer  Program  (UHDS) 
with  those  from  DATAPLOT  Software  Package  and  the  Measured  Values 


Calculated 

Va  lues 

Measured 

Va  lues 

Pine  No.  1 

Pine  No.  2 

Pipe 

UHDS  DATAPLOT 

UHDS  DATAPLOT 

No.  1 

No.  2 

Data  Set  No.  1: 

Heat  loss  rate , Btu/h-f t 

9,17 

7.11 

39.57 

41.56 

14.14 

43.89 

Horizontal  distance,  inch 

12.55 

12.61 

14.37 

14.42 

12.10 

13.92 

Vertical  depth,  inch 

10.05 

10.01 

10.13 

10.06 

10.66 

10.66 

Data  Set  No.  2: 

Heat  loss  rate,  Btu/h-ft 

13.00 

13.75 

68.61 

67.06 

11.71 

70.23 

Horizontal  distance,  inch 

12.07 

12.29 

13.88 

14.11 

12.10 

13.92 

Vertical  depth,  inch 

10.00 

9.93 

10.50 

10.45 

10.66 

10.66 
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Appendix  A,  The  Computer  Outputs  for  Sample  Calculations 


DATA  SET  NO.  1 


ITERATION 

NUMBER 

EPS 

PS  I 

RESIDUAL 
STD  DEV 

PARAMETER  ESTIMATES 
X(l)  TO  X(5) 

0 

7.22096 

14.000000 

12. 100000 
10.700000 

44.000000 
10.700000 

1 

1 . 00000 

. 68767 

5.81299 

11.224715 

12.439106 

10.414706 

41.630718 

10.628963 

2 

. 75000 

.28715 

5.61968 

9.458688 

12.615608 

9.304928 

40. 176298 
10.375861 

3 

1 . 12500 

. 24053 

5.49612 

9.243616 

12.613121 

10.644696 

40.056225 

10.256550 

4 

1.68750 

. 54323 

5.40505 

9. 174720 
12.607260 
9.940761 
40.019912 
10.217234 

5 

2.53125 

.41223 

5-39623 

9. 179133 
12.599882 

10. 130906 
40.016185 
10.201138 

6 

2-53125 

-51276 

5.39126 

9. 171416 
12.594186 
10.016180 
39.999598 

10. 189176 

7 

1 . 89844 

. 26886 

5.38884 

9. 163672 
12.585866 

10. 149976 
39.970063 

10. 172290 

8 

2.84766 

.47180 

5.38512 

9. 152807 
12.583094 
10.037032 
39.950680 

10. 166992 

9 

2.84766 

.70701 

5.38386 

9. 145402 

12.580213 
10.073288 
39.933222 
10. 162186 


9 


10 

11 

12 

13 

14 

15 

16 

17 

IB 

19 

20 


1 . 423B3 


57521 


2.40271 


2.40271 


2.70305 


2.70305 


2.02729 


3.04093 


3.04093 


1.52046 


2.5657B 


2.5657B 


5.3BIB7  9 

12, 

10 
39 
10 

.17232  5.38129  9 

12 

10 

39 

10 

.17407  5.38072  9 

12 

10 

39 

10 

.26952  5.38004  9 

12 

10 

39 

10 

.36568  5.37958  9 

12 

10 

39 

10 

.12593  5.37944  9 

12 

10 

39 

10 

.42682  5.37895  9 

12 

10 

39 

10 

.68200  5.37878  9 

12 

10 

39 

10 

.68249  5.37845  9 

12 

10 

39 

10 

,33706  5.37831  9 

12 

10 

39 

10 

.37661  5.37819  9 

12. 

10 

39, 

10 


115063 

572B72 

016444 

862524 

149941 

111916 

5702B1 

103933 

843351 

146056 

103438 

568392 

016042 

820004 

143633 

102192 

566727 

086757 

805844 

141446 

098956 

565406 

039324 

790031 

139923 

100236 

563155 

089440 

767473 

137136 

099056 

562354 

046445 

756064 

136357 

099242 

561498 

060547 

745748 

135552 

102387 

559048 

041647 

707701 

133321 

105828 

558151 

067994 

696178 

132491 

108172 

557429 

044659 

683856 

131969 


10 


21 

22 

23 

24 

25 

26 

27 

2B 

29 

3e 


1 . 92434 


05253 


2.88651 


2.88651 


1 . 44325 


2.43549 


2.43549 


2.73993 


2.73993 


3.08242 


3.08242 


5.37817  9 

12 
10 
39 
10 

.33563  5.37797  9 

12 

10 

39 

10 

.47197  5.37787  9 

12 

10 

39 

10 

.11607  5.37784  9 

12 

10 

39 

10 

.04983  5.37780  9 

12 

10 

39 

10 

.06086  5.37774  9 

12 

10 

39 

10 

.20914  5.37758  9 

12 

10 

39 

10 

.27602  5.37747  9 

12 

10 

39 

10 

.47820  5.37740  9 

12 

10 

39 

10 

.79349  5.37736  9 

12 

10 

39 

10 


116835 

556201 

077210 

665691 

130793 

119105 

555788 

044200 

656401 

130577 

122644 

555312 

060677 

648178 

130251 

135113 

554060 

032628 

615492 

129613 

142802 

553515 

077973 

606665 

129032 

147149 

553212 

030512 

595291 

128964 

153129 

552839 

070737 

588393 

128600 

157590 

552603 

043622 

580286 

128486 

161910 

552363 

058150 

574549 

128307 

166001 

552166 

052905 

568644 

128178 
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H-lBl 

i 

2 

3 

4 

5 

6 

7 

8 

9 

18 

11 

12 

13 

14 

15 

16 

17 

18 

19 

28 

21 

22 

23 

24 

25 

26 

27 

28 

29 

38 

31 

32 

33 

34 

35 

36 

37 


PARAMETERS 
X(  1) 

X(  2) 

X(  3) 

X(  4) 

X(  5) 


ESTIMATE 
9. 16680077 
12.55216611 
10.85290474 
39.56864357 
18. 12817812 


RESIDUAL  STANDARD  DEVIATION  = 5.37735825 

NUMBER  OF  RESIDUAL  DEGREES  OF  FREEDOM  * 


STANDARD 

DEVIATION 

8.72888001 

. 16923710 

2.54018502 

9,85351317 

,63176696 

43 


HORIZONTAL 
DISTANCE 
(IN.  ) 


4.088 

8.000 

10.000 

12.000 

14.000 

16.000 
18.000 
22.000 

4.000 

8.000 
10.000 
12.000 

14.000 

16.000 
18.000 
22.000 

4.000 

8.000 
10.000 
12.000 

14.000 

16.000 
18.000 
22.000 

4.000 

8.000 
10.000 
12.000 

14.000 

16.000 
18.000 
22.000 

4.000 

8.000 
10.000 
12.000 

14.000 

16.000 
18.000 
22.000 


VERTICAL 
DEPTH 
(IN.  ) 


1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

3.000 

3.000 

3.000 

3.000 

3.000 

3.000 

3.000 

3.000 

5.000 
5.000 
5.000 
5.000 
5.000 
5.000 
5.000 

5.000 

7.000 
7.000 
7.000 
7.000 
7.000 
7.000 
7.000 

7.000 

9.000 
9.000 
9.000 
9.000 
9.000 
9.000 
9.000 
9.000 


OBSERVED 
TEMP 
(DEG  F) 


85.80000 
87.20000 

87.90000 

88.50000 
90.60000 

90.40000 
89.70000 

88.00000 

89.70000 

92.30000 

94. 10000 

95.40000 

100.80000 

99.50000 

97.40000 

93.30000 

95.80000 

101.20000 

105.40000 

108.90000 

114.50000 

112.00000 

109.00000 

100.80000 

101.80000 

112.00000 

118.50000 

126. 10000 

132.50000 

129. 10000 

122.40000 
107.00000 

106.40000 

122.50000 

133.30000 

147.80000 

157.30000 

148.90000 

134.80000 

111. 90000 


PREDICTED 
TEMP 
(DEG  F) 


87,49305 
88.45967 
88.97559 
89.37932 
89.53824 
89.38827 
88.98901 
87.95742 
95. 19702 
98. 16691 
99.83375 
101.19498 
101.74833 
101.22743 
99.87854 
96.60349 
103.23637 
108,36755 
111. 56476 
114.45804 
115.74085 
114.54075 
111.65511 
105.59873 
110.28734 
117.61514 
122.95314 
128.86840 
132. 12347 
129. 10638 
123.09964 
113.53755 
115.70734 
124.70498 
132.37382 
144.59058 
156.35913 
144.95282 
132.51041 
119.58487 


RESIDUAL 
(DEG  F) 


-1,69305 
“1.25967 
“1.07559 
-.87932 
1.06176 
1.01173 
.71099 
. 04258 
-5.49702 
-5.86691 
-5.73375 
-5.79498 
-.94833 
-1.72743 
-2.47854 
-3,30349 
-7.43637 
-7. 16755 
-6. 16476 
-5.55804 
-1.24085 
-2.54075 
-2.65511 
-4.79873 
-8.48734 
-5.61514 
-4.45314 
-2.76840 
. 37653 
-.00638 
-.69964 
-6.53755 
-9,30734 
-2,20498 
.92618 
3.20942 
.94087 
3.94718 
2.28959 
-7.68487 
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41 

4. 008 

10.000 

108. 10000 

117.48310 

. -9.38310 

42 

8.000 

10.000 

127.40000 

126.84306 

. 55694 

43 

10.000 

10.000 

141.40000 

135.04996 

6.35004 

44 

12.000 

10.000 

163.70000 

151.16985 

12.53015 

45 

14.000 

10.000 

178.70000 

178.04528 

. 65472 

46 

16.000 

10.000 

161.70000 

150.25413 

1 1 . 44587 

47 

18.000 

10.000 

139.80000 

135. 18119 

4.61881 

48 

22.000 

10.000 

113.00000 

121.50745 

-8.50745 

PIPE  NO.  1 

PIPE  NO.  2 

HEAT  LOSS  RATE(Q) ,BTU/H-FT 

9. 1660 

39.5686 

HOR I ZONT AL  D I STANCE  < L) , I NCH 

12.5522 

14.3672 

VERTICAL  DEPTH (D) , INCH 

10.0529 

10. 1282 

13 


DATA  SET  NQ«  2 


ITERATION 

NUMBER 

EPS 

PS  I 

RESIDUAL 
STD  DEV 

PARAMETER  ESTIMATES 
X(l)  TO  X(5) 

0 

10.27628 

12.000000 

12. 100000 

10. 700000 

70. 000000 
10.700000 

1 

I . 00000 

. 08425 

10.20687 

1 1 . 894454 

12. 168653 
9.031536 
69.864176 
10,598553 

2 

.75000 

. 03475 

10. 18736 

12.352496 

12. 152267 
10.812442 
69.600809 
10.550360 

3 

1.12500 

o 05895 

10. 15070 

12.202898 

12. 172168 

9. 126657 
69.350817 
10.541380 

4 

1.12500 

.66719 

9.90287 

12.538536 

12, 148719 
10.248181 
69.350679 
10.531648 

5 

1 . 26563 

. 30277 

9.82288 

12.611932 

12. 145264 
9.841317 

69. 163855 
10.523090 

6 

1 . 26563 

. 35655 

9.76979 

12.764740 

12. 126046 
10.076785 
69.006133 
10.515387 

7 

2. 13574 

.48417 

9.73371 

12.819186 

12. 113955 
9.966080 
68.945657 
10.512808 

8 

3.20361 

. 46578 

9.72503 

12.852451 

12. 107520 
10.013358 
68.921964 
10.511583 

9 

7.20813 

. 24960 

9.72305 

12.859282 

12. 106218 
9.979526 
68.917197 
10.511344 


14 


10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


7.20B13 


79971 


9.72057 


12.866356 
12. 104903 
9.990806 
68.912568 
10.511102 


5.40610 


46249 


9.71916 


12.878060 
12. 102671 
10.008261 
68.903920 
10.510669 


12. 16372 


48388 


9.71778 


12.880325 
12. 102233 
9.991912 
68.902162 
10.510584 


12. 16372 


90221 


9.71728 


12.882628 
12. 101793 
9.995079 
68.900420 
10.510499 


6.08186 


99189 


9.71553 


12.891324 
12. 100083 
9.993929 
68.893189 
10.510159 


3.04093 


71147 


9.71095 


12.918940 

12.094011 

10.011531 

68.861410 

10.508773 


7.69735 


37750 


9.70883 


12.922895 

12.093095 

9.987494 

68.856139 

10.508560 


7.69735 


74186 


9.70752 


12.927037 

12.092173 

9.996239 

68.850985 

10.508344 


3.84868 


45828 


9.70584 


12.941556 

12.088684 

9.978789 

68.829463 

10.507472 


2.88651 


22830 


9.70277 


12.956534 

12.083693 

10.019330 

68.787850 

10.505857 


4.32976 


46995 


9.69722 


12.962753 

12.081517 

9.988662 

68.769102 

10.505136 


15 


21 


4.32976 


22 


23 


24 


25 


26 


27 


7.30647 


7.30647 


3.65324 


1 . 82662 


23.40711 


23.40711 


9.69599 


. 34749 


9.69434 


75312 


9.69321 


99299 


9.69042 


95221 


9.68187 


49086 


9.68161 


.99158 


9.68156 


12-969840 

12.079355 

10.009409 

68.751080 

10.504392 

12.972152 

12.078608 

9.988858 

68.744614 

10.504132 

12.974712 

12.077853 

9.996069 

68.738315 

10.503869 

12.983305 

12.075023 

9.995695 

68.712882 

10.502797 

12.998717 

12.066224 

10.001538 

68.612464 

10.498272 

12.998827 

12.066170 

9.996065 

68.611844 

10.498245 

12.998942 

12.066116 

9.996294 

68.611226 

10.498218 


PARAMETERS 
X(  1) 

X(  2) 

X(  3) 

X(  4) 

X(  5) 


ESTIMATE 

12.99894225 

12.06611595 

9.99629431 

68.61122629 

10.49821786 


RESIDUAL  STANDARD  DEVIATION  = 9-68155523 

NUMBER  OF  RESIDUAL  DEGREES  OF  FREEDOM  » 


STANDARD 

DEVIATION 

66.74597333 

1.19993271 

4.61363255 

68. 17682878 

. 46965895 

43 
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HORIZONTAL 

VERTICAL 

OBSERVED 

PREDICTED 

DISTANCE 

DEPTH 

TEMP 

TEMP 

RESIDUAL 

NUMBER  (IN.) 

(IN.  ) 

• (DEB  F) 

(DEG  F> 

(DEG  F) 

1 

4.888 

1.000 

86.40000 

90.25388 

-3.85388 

2 

8.888 

1.000 

91.10000 

91.84095 

-.74095 

3 

18.888 

1.000 

94.00000 

92.62902 

1.37098 

4 

12.880 

1.000 

96.00000 

93. 18403 

2.81597 

5 

14.000 

1.000 

95.00000 

93.32018 

1 . 67982 

6 

16.080 

1.000 

94.80000 

92.98568 

1.81432 

7 

18.808 

1.000 

93.70000 

92.30552 

1.39448 

8 

22.000 

1.000 

90. 10000 

90.67685 

-.57685 

9 

4.000 

3.000 

94.80000 

104. 13674 

-9.33674 

18 

8.000 

3.000 

102.00000 

109.05058 

-7.05058 

11 

10.000 

3.000 

105.80000 

111.61855 

-5.81855 

12 

12.000 

3.000 

120.00000 

113.49617 

6.50383 

13 

14.000 

3.000 

109.90000 

113.96602 

-4.06602 

14 

16.000 

3.000 

108.30000 

112.81577 

-4.51577 

IS 

18.000 

3.000 

105.30000 

110.54979 

-5.24979 

16 

22.000 

3.000 

98.40000 

105.42126 

-7.02126 

17 

4.000 

5.000 

100.30000 

115.50676 

-15.20676 

18 

8.000 

5.000 

114.30000 

124. 13741 

-9.83741 

19 

10.000 

5.000 

121.70000 

129.17252 

-7.47252 

28 

12.000 

5.000 

140.90000 

133.21193 

7.68807 

21 

14.000 

5.000 

129.80000 

134.27981 

-4.47981 

22 

16.000 

5.000 

126.60000 

131.69740 

-5.09740 

23 

18.000 

5.000 

121.50000 

127.01064 

-5.51064 

24 

22.000 

5.000 

108.50000 

117.67864 

-9. 17864 

25 

4.000 

7.000 

107.50000 

125.68172 

-18. 18172 

26 

8.000 

7.000 

128.50000 

138.32933 

-9.82933 

27 

10.000 

7.000 

140.60000 

147. 13762 

-6.53762 

28 

12.000 

7.000 

170.40000 

155.73634 

14.66366 

29 

14.000 

7.000 

156.00000 

158.36213 

-2.36213 

38 

16.000 

7.000 

151.20000 

152.26362 

-1.06362 

31 

18.000 

7.000 

139.80000 

143. 14326 

-3.34326 

32 

22.000 

7.000 

117.70000 

128.70848 

-11.00848 

33 

4.000 

9.000 

115.70000 

132.60822 

-16.90822 

34 

8.000 

9.000 

140.80000 

148.61248 

-7.81248 

35 

10.000 

9.000 

161.60000 

162.31705 

-.71705 

36 

10.000 

9.000 

188.20000 

162.31705 

25.88295 

37 

14.000 

9.000 

191.90000 

192.08542 

-. 18542 

38 

16.000 

9.000 

181.20000 

172.74064 

8.45936 

39 

18.000 

9.000 

158.00000 

155.68434 

2.31566 

48 

22.000 

9.000 

123.80000 

136.29532 

-12.49532 

41 

4.-000 

10.000 

118.60000 

134.97962 

-16.37962 

42 

8.000 

10.000 

146.20000 

151.83121 

-5.63121 

43 

10.000 

10.000 

173. 10000 

167.01626 

6.08374 

44 

12.000 

10.000 

212.20000 

206.44361 

5.75639 

45 

14.000 

10.000 

225.20000 

226.87275 

-1.67275 

46 

16.000 

10.000 

200.20000 

180.43853 

19.76147 

47 

18.000 

10.000 

165.70000 

159.65984 

6.04016 

48 

22.000 

10.000 

125.00000 

138.85673 

-13.85673 

PIPE  NO.  1 

PIPE  NO.  2 

HEAT  LOSS  RATE(Q) ,BTU/H-FT 

12.9989 

68.6112 

HORIZONTAL  DISTANCE (L) , INCH 

12.0661 

13.8811 

VERTICAL  DEPTH (D) , INCH 

9.9963 

10.4982 
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Appendix  B.  A Listing  of  the  Computer  Program 


PROGRAM  UHDS 

C MAIN  PROGRAM  FOR  UNCONSTRAINED,  UNWEIGHTED  NONLINEAR  LEAST  SQUARES 

C FITTING  USING  THE  LEVENBERG/MARQUARDT/MORRI SON  ALGORITHM  WITH 

C ANALYTICAL  DERIVATIVES.  THESE  DIMENSIONS  ALLOW  UP  TO  100 

C EQUATIONS  IN  10  UNKNOWNS  WITH  UP  TO  5 INDEPENDENT  VARIABLES, 

C SUBROUTINES  CALLED  INCLUDE  LMMNLF  AND  FUNVAL. 

C INPUT  FILE  IS  DATAFL 

IMPLICIT  REAL*8  (A-F,S“Y) 

DIMENSION  X(10),YY(100),XX(100,5),F(100),A(100,10) 

COMMON  YY,XX,NSET 

C** 

IER*=2 
ITS-=50 
TOL“l ,D-5 
EPS*l.D-8 
EXPEND«1 .5 
DECR-0.5 

C**  NUMBERS  OF  INDEPENDENT  VARIABLES  AND  PARAMETERS  IN  THE 
C**  NONLINEAR  EQUATIONS 

NIV«5 
NP*5 

OPEN(8,FILE*^DATAFL' ) 

' OPEN(6,FILE-'PRN') 

C READ  IN  THE  NUMBER  OF  DATA  SETS 

READ(8, 5, END-100, ERR“150)  NSETS 
5 FORMAT  (15) 

NSET-1 

C READ  IN  THE  NUMBER  OF  DATA  POINTS,  SOIL  THERMAL  CONDUCTIVITY,  AND 

C THE  SEPARATION  DISTANCE  BETWEEN  THE  CENTERS  OF  THE  PIPE 

10  READ(8,15,ERR-150)  N,AK,DS 

15  FORMAT  (I5,2F10.5) 

WRITE(6,18)  NSET 

18  F0RMAT(1H1,1X,'DATA  SET  NO. ',13/) 

C READ  IN  THE  MEASURED  EARTH  TEMP,  THE  HORIZONTAL  DISTANCE,  THE 

C VERTICAL  DEPTH.  AND  THE  UNDISTURBED  EARTH  TEMP 

DO  30  I-1,N 

READ(8,20,ERR-150)  YY(l),  (XX (I , J ) , J-1 , 3 ) 

20  FORMAT  (F5.1,2(lX,F4.1),lX,F5.1) 

30  CONTINUE 

C COMPUTE  RH0-1/(4*PAI*K) ,AND  THE  DISTANCE  BETWEEN  THE  PIPE  CENTERS 

DO  40  I-1,N 

XX(I ,4)-l./(4.*3.14159*AK) 

XX(I ,5)-DS 
40  CONTINUE 

C READING  IN  THE  INITIAL  PARAMETER  ESTIMATES 

C THE  PARAMETERS  ARE  Q1 , B 1 ( H . DI STANCE ) ,A1 ( PIP E DEPTH), Q2,  AND  A2 

DO  60  I«1,NP 
READ(8,50,ERR-150)  X(I) 

50  FORMAT  (F10.5), 

60  CONTINUE 

CALL  LMMNLF  (X , F , A , SUMSQ , N , NP , TOL , EXP END , DECR , I TS , I ER ) 

IF  (lER  .EQ.  2)  WRITE(6,80) 

80  FORMAT  (IX,'  MAXIMUM  NUMBERS  OF  ITERATIONS  EXCEEDED  ') 

C PRINT  THE  HEAT  LOSS  RATES  FROM  THE  UNDERGROUND  PIPES  AND  THEIR 

C LOCATIONS 
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WRITE(6,90) 

90  FORMAT(//36X,'  PIPE  NO.  1 " , 6X , ' PIPE  NO.  2'/) 

DHL2«X(2)  + DS 

WRITE(6 ,95)  X(l) ,X(4) ,X ( 2 ) , DHL2 ,X ( 3 ) , X ( 5 ) 

95  FORMAT(2X,'HEAT  LOSS  RATE ( Q ) , BTU /H-FT" , 2 ( 8X , F 1 0 . 4 ) / 2X , ' HORI ZONTAL 
& DISTANCE(L) .INCH' , 6X , F 1 0 . 4 , 8X , F 1 0 . 4/ 2X , ' VERTI CAL  DEPTH(D) ,INCH' , 
&12X,F10.4,8X,F10.4) 

NSET*NSET+1 

IF  (NSET  .LE.  NSETS)  GOTO  10 
100  WRITE(*,120) 

120  FORMAT ( IX, 'DONE' ) 

GOTO  300 

150  WRITE(*,160) 

160  FORMATC IX, 'THERE  IS  AN  ERROR  IN  INPUT  DATA') 

300  CLOSE(8) 

STOP 

END 


19 


o o 


SUBROUTINE  FUNVAL  (A  ,F ,X , SUMSQ ,IFL ,N) 

THIS  SUBROUTINE  IS  USED  WITH  SUBROUTINE  LMMNLF  TO  EVALUATE  TEE 
FUNCTION  G AND  ITS  DERIVATIVES. 

IMPLICIT  REAL*8  (A-G,R-Y) 

DIMENSION  X(10),YY(100),XX(100,5),F(100),A(100,10) 

REALMS  NUMl , NUM2 , NUM3 , NUM4 , NUM5 , NUM6 
COMMON  YY,XX,NSET 
SDMSQ“0.D0 
DO  10  I-1,N 

NUMl* (XX (I ,1)-X(2))**2+(XX(I ,2)+X(3))**2 
DENI* (XX (I ,1)-X(2) )**2+(XX(l ,2)-X(3) )**2 
NUM2-(XX(I ,1)-X(2)-XX(I , 5 ) ) **2+ (XX (I , 2 ) +X ( 5 ) ) **2 
DEN2«(XX(I ,1)-X(2)-XX(I ,5))**2+(XX(l ,2)-X(5) )**2 
NUM3*XX(I ,2)*(XX(I ,1)-X(2)) 

DEN3*NUM1*DEN1 

NUM4*XX(I »2)*(XX(I ,1)-X(2)-XX(I ,5)) 

DEN4«NUM2*DEN2 

NUM5*XX(I ,2)*( (XX(I , 1)-X(2) )**2+(XX(I , 2 )**2-X ( 3 )**2 ) ) 

NUM6*XX(I ,2)*( (XX(I ,1)“X(2)-XX(I , 5 ) )**2+ (XX ( I , 2 )**2-X ( 5 )**2 ) ) 

C CALCULATE  THE  VALUE  OF  FUNCTION  G 

G*XX(I ,4)*(X(1)*DL0G(NUM1/DEN1)+X(4)*DL0G(NUM2/DEN2) )+XX(l ,3) 
RESID*YY(I)~G 
SUMSQ*SUMSQ+RESID*RESID 
IF  (IFL  .EQ.  2)  GOTO  10 
C SET  VALUES  FOR  I-TH  ROW  OF  GRADIENT  G 

A(I ,1)*-XX(I ,4)*DL0G(NUM1/DEN1) 

A(I  ,2)— 8.*XX(I  ,4)*‘(X(1)*X(3)*NUM3/DEN3+X(4)*X(5)*NUM4/DEN4) 
A(I ,3)--4.*XX(I ,4)*X(1)*NUM5/DEN3 
A(I ,4)*-XX(l ,4)*DL0G(NUM2/DEN2) 

A(I ,5)--4.*XX(I ,4)*X(4)*NUM6/DEN4 
F(I )*RESID 
10  CONTINUE 

RETURN 
END 
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SUBROUTINE  LMMNLF  ( X , F , A , SUMSQ , N , NP , TOL , EXPND , DECK , I TS , I ER) 
IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  B(10,10),DA(10),DU(10),D(10),C(10),DX(10),Y(10) 

DIMENSION  X(10) ,YY( 100) ,XX( 100,5 ) ,F( 100) ,A( 100, 10) 

COMMON  YY,XX,NSET 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


THIS  SUBROUTINE  IS  BASED  ON  LEVENBURG,  MARQUARDT,  MORRISON 
ALGORITHM  (SEE  OSBORNE  'NONLINEAR  LEAST  SQUARES  - THE  LEVENBERG 
ALGORITHM  REVISITED',  J.  AUSTRAL.  MATH.  SOC.  19  (SERIES  B)  (1976), 


PP.  343-357)  AND 
IN  THE  NONLINEAR 


IS  MODIFIED 
FUNCTION. 


FOR  ONE  OR  MORE  INDEPENDENT  VARIABLES 


VARIABLES 

X(l) 


A(N,NP) 


Vector 

Input 

Output 

Matrix 
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F(l) 


of  parameters  less  than  or  equal  to 
: Contains  estimate  of  solution 
: Contains  solution  vector 

containg  the  first  partial  derivatives  of  the  functi 
with  respect  to  each  of  the  parameters. 

Output  : Contains  Upper  Triangular  Factor  in  orthogonal 
factorization  of  GRAD  F 

Storage  for  F vector  of  terms  in  sum  of  squares 


c 

SUMSQ 

Output 

: Contains 

final  residual  sum  of  squares 

c 

N 

Input 

: Dimension 

of  F 

c 

NP 

Input 

: Dimension 

of  X 

c 

TOL 

Input 

: Tolerence 

on  Calculation 

c 

EXPND 

Input 

: Factor  by 

which  EPS  increased  if  test  on 

sum  0 f 

c 

squares 

f *a  i 1 8 

c 

DECR 

Input 

: Factor  by 

which  EPS  decreased  if  test  on 

sum  of 

c 

squares 

succeeds  on 

first  attempt 

ITS 


lER 


Input 
Output  : 
Input  : ■ 
-1 

Output  : *1 

-2 

-3 

-4 

If  lER  *2,3 


Max  number  of  iterations 
Actual  number  of  iterations 
*0  No  Printing 

Print  Diagnostic  Information 
Successful  Termination 
Max  ITS  Exceeded 
EPS  exceeds  1.D6 


or  4 


Reached 

Tol 

too 

sma  1 1 

rors  in 

grad 

ien  t 

as  a sea 

le  which 

is 

an  norm 

of  A 

by 

a 

0 set  va 

lues 

of 

SUMSQ, 

calculat ion 

*500+1  I'th  column  of  A 
small  compared  to  Eucli< 

Factor  less  than  1.D6 
User  supplied  subroutine  FUNVAL  required 
F,  A.  Declaration  must  be 

SUBROUTINE  FUNVAL  ( A , F ,X , SUMSQ , IFL , N ) 

If  IFL*1  sets  all  values 

If  IFL*2  sets  SUMSQ  only;  must  not  alter  A or  F 
Diagonostic  information  contains  in  an  output  file:  DIAGON.DTA 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

NRDF*N-NP 

IPRINT*IER 

IF  (IPRINT .EQ. 0)  GO  TO  41 
IF(NSET.GT.1)G0  TO  300 
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,STATUS*'NEW' ) 


OPEN  (3,FILE»='DIAGON.DTA  ' 

WRITE  (3,102) 

41  MAXITS-ITS 

300  WRITE(3,190)  NSET 

WRITE(6 ,200) 

ITSO-0 

CALL  FUNVAL(A,F,X,SSF, 1 ,N) 

SDRES-DSQRT(SSF/NRDF ) 

WRITE(6 ,201)ITS0 ,SDRES ,X(1) 

DO  210  1-2, NP 
WRITE(6,202)  X(I) 

210  CONTINUE 
ITS-O 

40  ITS-ITS+1 

NITS-0 

CALL  FUNVAL(A,F,X,SSF,1 ,N) 

C COMPUTE  ESTIMATE  OF  RESIDUAL  STANDARD  DEVIATION 

CCCCCCCCCCCCCCCCCCCCCCCC 
C SCALE  GRAD  F C 

CCCCCCCCCCCCCCCCCCCCCCCC 
W-O.DO 
DO  1 1*1, NP 
S*0,D0 
DO  2 J*1  ,N 

2 S*S+A(J,I)**2 
W-W+S 

1 D(I)*DSQRT(S) 

W*DSQRT(W) 

DO  46  1*1, NP 

IF  (D(l)/W,LT.l.D-6)  GO  TO  47 
S*1.0/D(I) 

DO  3 J»1,N 

3 A(J,I)-A(J,I)*S 

46  CONTINUE 
GO  TO  48 

47  IER-500+I 

IF  (IPRINT.EQ.O)  GO  TO  49 

WRITE(3,104)  I 

WRITE(3 ,105)  (D(I),I*1,NP) 

49  GO  TO  45 

48  IF  (ITS.EQ.l)  EPS-1.0 

IF  (IPRINT.EQ.O)  GO  TO  42 
WRITE(3,100)  ITS,EPS,SSF 

ccccccccccccccccccccccccccccccccccccccccccccccccccc 

C HOUSEHOLDER  TRANSFORMATION  OF  GRAD  F,F  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCQCCCCCCCCCCCCCCC 
C VECTOR  DA  CONTAINS  DIAGONAL  ELEMENTS  OF  UPPER 

C TRIANGULAR  MATRIX  A. 

42  DO  4 1*1, NP 
S*0oD0 

DO  5 J*I,N 
5 S*S+A(J,I)**2 

S*DSQRT(S) 

IF  (A(I  ,1)  .GT.0.0)  S — S 
DA(I)-S 
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A(I  ,I)»A(I ,I)-S 
IF  (I.EQ.NP)  GO  TO  6 
IPl-I+1 
DO  7 K«IP1,NP 
S“0.D0 
DO  8 J-I,N 
S-S+A(J,I)*A(J,K) 

S — S/(DA(I)*A(I  ,1)) 

DO  9 J-I ,N 
9 A( J,K)«A( J ,K)-S*A( J ,I ) 

7 CONTINUE 

6 S-O.DO 

DO  20  J»I,N 

20  S»S+A(J,I)*F(J) 

S— S/(DA(I)*A(I,I)) 

DO  21  J-I.N 

21  F(J)-F(J)-S*A(J,I) 

4 CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C COMPUTE  SUM  OF  SQUARES  OF  RESIDUALS  C 

ccccccccccccccccccccccccccccccccccccccccccccccc 

, NP1*NP+1 
SSR>=0.D0 
DO  22  I«NP1,N 

22  SSR»SSR+F(I)**2 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C FACTOR  EPS  APENDAGE , TRANSFORM  RES  UPPER  TRIANGLE  OF 

C TRANSFORMED  MATRIX  STORED  IN  UPPER  TRIANGLE  OF  B. 

C FILL  IN  B STORED  COLUMNWISE  IN  ROWS  IN  LOWER  TRIANGLE  OF  B. 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

19  DO  30  1*1, NP 

DO  31  J*1,NP 
31  B(I,J)-0.D0 

C(I)«0.D0 
30  B(I,I)*EPS 

DO  10  1*1, NP 
S*DA(I)**2 
IP1*I+1 
IL1*I-1 
DO  12  J*1,I 
12  S*S+B(I ,J)**2 

S*DSQRT(S) 

IF  (DA(I)  .GT.O.DO)  S — S 

DU(I)*S 

W*DA(I)-S 

IF  (I.EQ.NP)  GO  TO  18 
DO  13  K*IP1 ,NP 
S*A(I ,K)*W 

IF  (I.EQ.I)  GO  TO  11 
DO  14  J*1,IL1 

14  S*S+B(I  ,J)*B(K, J) 

11  S*-S/(DU(I)*W) 

B(I ,K)*A(I ,K)-S*W 
DO  15  J*1,I 

15  B(K, J)*B(K, J)-S*B(I ,J) 
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ODD 


13  CONTINUE 

18  S-F(I)*W 

DO  16  J*1,I 

16  S*S+B(I ,J)*C( J) 

S«-S/(DD(I )^W) 

DX(I)-F(I)-S*W 
DO  17  J-1,I 

17  C( J)-C(J)-S*B(I , J) 

10  CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C BACK  SUBSTITUTION  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DX  ( NP  ) *DX  ( NP  ) / DU  ( NP  ) 

DO  25  1*2, NP 
K*NP-I+1 
KPl-K+1 
S-O.DO 

DO  26  J*KP1,NP 
26  S*S+B(K,J)*DX( J) 

25  DX(K)*(DX(K)-S)/DU(K) 

SSS-SSR 

^ DO  32  1*1, NP 
SSS*SSS+C(I )**2 
DX(I)-DX(I)/D(I) 

32  Y(I)*X(I)-DX(I) 

NITS*NITS+1 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C CHECK  CONVERGENCE  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
IER«4 

IF  (SSS.GE.SSF)  GO  TO  45 
IER*1 

CALL  FUNVAL(A,F,Y,SSN,2,N) 

S-.5D0*(SSF-SSN)/(SSF-SSS) 

IF  (IPRINT.EQ.O)  GO  TO  43 

43  IF  (S.GE.l,D-4)  GO  TO  28 
EPS*EXPND*EPS 

IEE*3 

IF  (EPS.GT.1.D6)  GO  TO  45 
GO  TO  19 

28  SDRES-DSQRT(SSN/NRDF) 

DO  29  1*1, NP 

29  X(I)*Y(I) 

IF  (IPRINT.EQ.O)  GO  TO  44 
WRITE(6 ,203)  ITS , EPS , S , SDRES ,X(1) 

DO  211  1*2, NP 
WRITE(6,202)  X(I) 

211  CONTINUE 

C CHECK  FOR  CONVERGENCE  OF  SUM  OF  SQUARES  OF  RESIDUALS. 

44  IF  ((DSQRT(SSF)-DSQRT(SSS))/(1.D0  + DSQRT ( S SF ) ) . GE . TOL ) GO  TO  35 

45  SUMSQ-SSN 

DO  33  1*1, NP 
A(I ,I)*DA(I) 

S-D(I) 

DO  34  J*1,I 
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34  A(J,I)-A(J,I )*S 

33  CONTINUE 

C PRINT  ESTIMATES  OF  PARAMETERS  AND  THEIR  STANDARD  DEVIATIONS 

WRITE(6 ,204) 

DO  270  K«1,NP 

51- O.DO 
D(K)-1/A(K,K) 

S1*=S1+D(K)**2 

KPl-K+1 

DO  260  I-KPl ,NP 

52- 0.D0 
ILl-I-1 

DO  250  J*K,IL1 
250  S2-S2+A(J,I)*D(J) 

D(I)— S2/A(I  ,1) 

S1«S1+D(I)**2 
260  CONTINUE 

S1-SDRES*DSQRT(S1 ) 

WRITE(6,265)  K,X(K) ,S1 
270  CONTINUE 

C PRINT  RESIDUAL  STANDARD  DEVIATION  AND  DEGREES  OF  FREEDOM 

' WRITE(6,266)  SDRES,NRDF 

C PRINT  OVSERVATIONS , PREDICTED  VALUES  AND  RESIDUALS 

WRITE(6,275) 

CALL  FUNVAL(A,F,X,SUMSQ,1,N) 

DO  280  I-1,N 
PRED-YY(I)-F(I) 

WRITE(6 ,276)  I ,XX ( I , 1 ) ,XX( I , 2 ) , YY ( I ) ,PRED,F(I) 

280  CONTINUE 
RETURN 

35  IER-2 

IF  (ITS.GE.MAXITS)  GO  TO  45 
IF  (NITS.EQ.l)  EPS«EPS*DECR 
GO  TO  40 

100  FORMAT  ('  ITS*^,I3,'  EPS*' ,F14. 6 SUMSQ*' ,F14 . 6 ) 

102  FORMAT  ('1  NONLINEAR  LEAST  SQUARES  FIT  BY  LEVENBERG  ALGORITHM') 

104  FORMAT  ('SCALING  ERROR  NO.  OF  COLUMN  *',I3) 

105  FORMAT  (4('  D ( ' ,I 2 , ' ) -' ,F 14 . 6 ) ) 

190  FORMAT (IX, 'DATA  SET  NO. ',13/) 

200  FORMAT  ( 2X ,' I TERATI ON' , 27X ,' RES IDUAL' , 5X ,' PARAMETER  ESTIMATES '/ 3X , 
&'NUMBER' ,7X,'EPS' ,9X,'PSI' ,7X,'STD  DEV' , 9X , 'X ( 1 ) TO  X(5)'/) 

201  F0RMAT(4X,I3 ,29X,F10.5,4X,F15.6) 

202  , FORMAT(50X,F15.6) 

203  F0RMAT(/4X,I3 ,5X,F10.5,2X,F10.5,2X,F10.5,4X,F15.6) 

204  FORMAT(//57X,' STANDARD' /I OX, 'PARAMETERS' , 13X ,' ESTIMATE' , 16X, 
S'DEVIATION'/) 

265  FORMAT (12X,'X(' ,I2,')',12X,F14.8,12X,F12.8/) 

266  FORMATdOX, 'RESIDUAL  STANDARD  DEVIATION  * ' , F 1 4 . 8 / / 1 OX  , ' NUMBER  OF 
&RESIDUAL  DEGREES  OF  FREEDOM  * ',18//) 

275  FORMATdOX, 'HORIZONTAL' , 3X , ' VERTI CAL' , 5X ,' OBSERVED' , 5X , ' PREDI CTED 
&' /I IX, 'DISTANCE' ,5X, 'DEPTH' ,9X,'TEMP' ,9X,'TEMP' , 8X ,' RESIDUAL' / 2X 
&, 'NUMBER', 4X, '(IN. )',7X, '(IN. )',7X, '(DEG  F)',6X,'(DEG  F)',7X, 
&'(DEG  F)'//) 

276  FORMAT (3X,I3 ,2(4X,F8.3) ,3(4X, FI 0.5) ) 

END 
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